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Abstract 

A new integration scheme, combining the stability and the precision of usual 
pseudo-spectral codes with the locality of finite differences methods, is in- 
troduced. It turns out to be particularly suitable for the study of front and 
disturbance propagation in extended systems. An application to the complex 
Ginzburg-Landau equation shows the higher precision of this method with 
respect to spectral ones. 
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I. INTRODUCTION 



The integration of partial differential equations (PDE's) is a relevant topic by itself, but 
efficient and accurate algorithms are particularly important for extended systems exhibit- 
ing spatio-temporal chaos. Usually, PDE's are integrated according to one of two general 
schemes: finite difference methods or time splitting pseudo-spectral codes [f]. In the present 
contribution we will introduce a new algorithm that will combine the accuracy and stability 
of the spectral method with the locality property of finite differences schemes. 

The present method can be applied to any PDE, be it Hamiltonian or dissipative. Eor 
most problems, it should be comparable in efficiency to time splitting pseudo-spectral codes, 
but there is one class of problems where it seems superior to any previously known method. 

This is the propagation of disturbances or fronts into unstable states. The classic example 
is the propagation of fronts in the KPP (or Eisher-) equation [2]. Similar are propagation 
of solidification fronts in directional solidification and viscous fingering [3,4]. Another phe- 
nomenon which is closely related mathematically though it seems at first rather different, 
is propagation of disturbances in spatially extended chaotic systems [5,6]. In these cases, 
spectral methods cannot be used because of the instability of the state into which the front 
or the disturbance propagates. Due to the non-locality of the Eourier transform, it is in 
general impossible to keep the state completely unperturbed in front of the perturbation. 
Instead, there will be small perturbations in regions which should not have been reached 
yet by the front (at least due to round-off errors), and these perturbations will grow expo- 
nentially due to the instability of the state, rendering any precise measurement of the front 
properties impossible. 

Eor sake of brevity, we will discuss in this article only an application to propagation of 
disturbances in the one-dimensional complex Ginzburg-Landau equation (CGLE). We will 
consider the I-d CGLE because it represents a paradigmatic example for a large class of 
spatially extended systems [3] exhibiting, for different parameter intervals, behavior that 
ranges from spatio-temporal chaos [7] to stable and intermittent regimes [8]. Moreover, 
also the nonlinear Schrodinger equation can be recovered as a particular limit of the CGLE. 
Propagation turns out to play a crucial role in the understanding of the dynamics of spatially 
extended chaotic systems [4-6,9]. It is particularly important in systems with so-called 
"chaotic supertransients" [10], where it is responsible for the chaos observed in infinite 
systems. 

The paper is organized as follows. In the next section the integration scheme will be 
described in detail together with a specific implementation of the algorithm for the CGLE. 
Some accuracy tests are also reported in Sec. II. Section III is devoted to a comparison 
between propagation speeds measured with our algorithm, and speeds obtained from the 
Lyapunov analysis of [6]. Einally, Sec. IV summarizes the main results of the work, along 
with a few concluding remarks. 

II. THE INTEGRATION ALGORITHM 

In the present article we will limit ourselves to PDE's of the form 
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'^^^ = CAix,t) = ir+V) Aix,t) (1) 

where A[x^t) is a 1 component field, and x G [0, is 1-dimensional. We assume that 
the operator C can be split into two parts, a linear operator T containing all the spatial 
derivatives and a strictly local operator V which is in general non-linear. Both T and V 
are assumed to be translation invariant. Generalizations to multi-component fields and to 
higher dimensions are straightforward. 

Formally the solution of Eq. (1) can be expressed as 

A{x,t) = e'^A{x,0) . (2) 

We assume that T and V do not commute with each other. Thus, one cannot factorize 
the exponential, but should use the Baker-Campbell-Hausdorff (BCH) formula to express 
gt(T+V) ^ product of terms e*'''^ and e*^'^. An appropriate choice of the coefficients di 
and Ci can reduce noticeably the errors committed in neglecting higher order commutators 
in the BCH expression [11-13]. For sake of simplicity, we will consider the lowest order 
approximation [14] corresponding to the Trotter formula 

This belongs to the class of the so-called time splitting algorithms [1] and it is commonly 
referred as "leap-frog method". In practice, when performing this repeatedly, it is advan- 
tageous to combine the two half steps in all but the first and last iterations. Fxcept for 
the latter, each iteration bringing us from t = nr to t' = [n -\- l)r corresponds then to two 
successive integration steps, 

An+i/2ix) = e^^Anix) ; An+lix) = e^'^A„+i/2(x) ; (4) 

The first step, containing no spatial derivatives, can easily be done in coordinate space 
where it consists just of a multiplication or of the solution of a simple ODF. The second 
step would be non-local in coordinate space, and is thus usually done in momentum space. 
Thus, the field A is first Fourier transformed so that the integration takes the form of a 
simple multiplication, 

= e^^(^)i„(p) , (5) 

and an inverse Fourier transform brings the field back into coordinate representation. 

The novelty of the our method, compared to the usual time splitting pseudo-spectral 
scheme described above, consists just in performing the second step also in configuration 
space, avoiding thereby the two Fourier transforms. The rest of the integration algorithm is 
left unchanged. If we neglect spatial discretization, we would have 



A„+i(x)= / </>(e)A„+i/2(x-0 (6) 

^ — oo 

with 
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m = 7^ [dp e'-^«+^^(^) . (7) 

ZTT J 

The crucial observation is that will be strongly centered at if the step r is small. 
Thus, when using a spatial grid, the integration can be replaced by a sum which involves 
only very few terms. Let us assume that we can neglect outside an interval [— 5*, 5*], and 
that we work with a grid size Ax. Then each application of eq.(6) involves Nc = {2S ■ Ax 
operations (the values of (j)[iAx) are of course computed only once during the initialization 
of the routine). This is to be compared to the 0{2ln[L/ Ax]) operations for the two FFT's 
in the spectral method. 

To be specific, let us now consider the CGLE 

At = {l + ici)A,, + A-{l-ics)\A\^A (8) 

where the field A[x^t) is complex, whereas the parameters Ci and C3 assume real positive 
values. The analysis reported in this paper will be limited to a parameter region where it is 
known that Eq. (8) will show a chaotic behaviour [7] (namely, Ci = 3.5 and 0.6 < C3 < 0.9). 

The integration involving the operator V can be performed rewriting the field as A[x^ t) = 
p[x^t)e''^^^'^'^ and solving the ordinary differential equations 

— — — = cs p {x,t) (9) 

and 

^^-^ = 2[p^{x,t)-p\x,t)]. (10) 

The solution of eq.(lO) is given by 

p{x,T) = [e-'^{l/p{x,OY-l) + l]-'/\ (11) 

Inserting this into eq.(9) gives 

V'(x, r) = i^ix, 0) + C3 {r + ln[p(x, 0)/p(x, r)]} . (12) 

Since T[p) is quadratic in p, the kernel for the second integration step reads simply 

(j){x) = y^e"^''"' [cos(Ax2) - I sin(Ax2)] (13) 

where fi = fir — ifii = (1 — zci)/[4r(l + c\)\. In fig.l the real and complex part of (t){x) are 
shown for a typical choice of parameters. It is clear that the shape of (t){x) will strongly 
depend on the chosen time step and on the parameter Ci in the CGLE. In order to achieve 
sufficient accuracy, the shape and the details of (t){x) should be well resolved within the 
chosen spatial resolution dx. In fig. 2 we show average errors accumulated during a fixed 
total integration time T = 0.2 for several values of r and of the number Nc of terms in the 
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convolution. In these simulations, the spatial resolution and the parameters of the CGLE 
were kept fixed at Ax = 0.098, Ci = 3.5 and C3 = 0.9. 

We see in particular that the errors saturate as Nc increases, and we verified that the 
limits of the errors for large Nc coincide exactly with the errors of the spectral method, 
as it should be the case. Assume that we want to tolerate a total error which is, say, a 
factor 2 larger than that of the spectral method. Then we can read off from fig. 2 that the 
necessary values of are ^ 57 (for r = 0.01), 70 (for r = 0.02), 90 (for r = 0.05), and 110 
(for T = 0.1). On the other hand, using the FFT routine C06ECF of the NAG library and 
grid sizes N which are powers of 2, we found that our algorithm (whose speed is essentially 
proportional to Nc) needed more CPU time by a factor 0.15A^c/ In A^, for A^ > 2048. Thus 
it seems that our algorithm with r = 0.01 and Nc = 57 has about the same speed and 
precision as the spectral code with r = 0.014, for grids with A^ = 2^^. Thus both algorithms 
are comparable for typical problems. 



III. FRONT PROPAGATION SPEED 



In this Section we will consider only the propagation of disturbances into a chaotic one 
dimensional system. As we said, essentially the same problems arise for front propagation 
into unstable states. Also, the extention to multi-dimensional system with steady or chaotic 
states would be straightforward. 

In order to measure the propagation of a disturbance, we consider two realizations of 
the same field, say A[x^t) and B{x,t), which differ only in a finite region of space near the 
origin. If these states are unstable, the difference AA(a;,t) = \A{x,t) — B{x,t)\ will spread 
with a limiting velocity vp. This velocity can be measured by recording the distance R{t) 
of the front. This is the furthest point where AA(a;,t) is not exactly equal to zero, 

VF = lim (14) 

t^oo t 

This definition seems difficult to implement in practice, because we can only study finite 
systems. If we use periodic boundary conditions, the two opposing fronts (with velocities 
Vp and —Vp) will collide with each other after L/2vf time units. With open boundary 
conditions, the situation is not much better. In order to avoid this problem, we proceeded 
as follows. Assume that we have discretized space and time, and we know that the front 
cannot proceed by more than m sites per time step (in our method, m = (2A^c — l)/2)- Then 
we take configuration A as reference configuration, and replace in each step B[x^t) by 

B{x,t) ^ A{x,t) Vx G ]i?(t),i?(t) + m] . (15) 

By "cleaning" the solution B in the region ahead of the front, we make sure that the front 
can propagate unperturbed, while the effect of this cleaning on the trailing end of the front 
should be negligible for sufficiently large systems. 

Let us now assume that we want to use a spectral method. Due to the non-local nature 
of the Fourier transform, the difference AA[x^ t) will not remain identically zero ahead of the 
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front. This means that the perturbation will spread instantaneously, albeit with very small 
amplitude (typically of the order of the machine precision). Since the states are unstable, 
this tiny error will increase exponentially, preventing thereby any accurate measurement of 

Speeds measured with our algorithm are given in table I, for Ci = 3.5 and several values 
of C3 (third column). Together with the speed, we measured also the average profile of the 
front of AA[x^t) in a co-moving coordinate system. Of course, since both A and B are 
chaotic, also AA will fiuctuate. But we expect [5,6] that log AA(a;,t) will decay in average 
exponentially in space, 

AA(x,t) ~ e"^^^ . (16) 

Estimates for the exponent fip are given in the second column of table I. 

Recently, it has beeen shown for coupled map lattices [15] that fip and vp are related 
[6]. More precisely, let us define the specific Lyapunov exponent A(/i) [16] as the average 
growth in time, at fixed position x, of an infinitesimal perturbation 5A[x^t) which satisfies 
5A[x^t) ~ e~^^ throughout the entire system (this can be achieved by special boundary 
conditions, SA{L,t) = e~^^SA{0,t)). Furthermore, we define 

V(,0 = ^ . (17) 

Then we have [6] 

VF = VM . (18) 

In addition, we know that V{fJ,) has a unique minimum for chaotic systems, so that [6] 

vp ^ vl = minV{ji) . (19) 

The speed vl is called the linear velocity^ since it is obtained from a linear stability analysis. 

The values of vl and of fip (the value of /i where V{fJ,) has its minimum) are given 
in the last two columns of table I. We see that they agree within errors with vp and /ij?, 
respectively. Some discrepancies are observed between the spatial profiles fip and fip (in 
particular, for C3 = 0.9), due to the difficulties to perform a direct measurement of the shape 
of the propagating front. 

We see thus that disturbances are "pulled" [17] for the present parameter values of the 
CGLE, except maybe for the smallest values of C3. That is, it is the very edge of the 
perturbation front (where AA is infinitesimal and described by linear stability analysis) 
which "pulls" it. In contrast, there are cases possible where the front is "pushed" by regions 
where the perturbation is finite. We will present in a forthcoming paper [18] evidence for 
pushed error propagation in the CGLE at different parameter values. 

IV. CONCLUSIONS 

In the present paper we have introduced a new technique to integrate PDE's. It consists 
of a time-splitting procedure (as, e.g., the well known leap frog) where the integration of 
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the operator which is non-local in configuration space is performed by convolution with its 
kernel. 

This algorithm is particularly useful for propagation phenomena into unstable (and 
chaotic, in particular) states which cannot be treated by pseudospectral methods. We ap- 
plied it to the propagation of disturbances in the CLGE, but we suggest that it will be useful 
also in many similar phenomena. 

The speed of disturbance propagation is a relevant indicator to characterize spatially 
extended systems, being able not only to account for local chaoticity (associated to a maximal 
positive Lyapunov exponent) but also for the information spreading in the system. In 
the present paper, the situation has been considered where the propagation is essentially 
ruled by pulling (i.e., by a linear mechanism). This seems not to be always the case. The 
possibility of having different propagation mechanisms should be closely related to various 
phases observed for the CLGE (in particular, to the two different chaotic regimes found 
for this system, namely phase and defect turbulence) [7], and should help to explain the 
dynamics in this model. Work along these lines is in progress and will be reported elsewhere 
[18]. 

ACKNOWLEDGMENTS 

We want to thank A. Politi for useful suggestions and discussions, as well as D.A. Egolf 
and H.S. Greenside for providing us some unpublished data concerning the CGLE. One of 
us (A.T.) gratefully acknowledges the European Economic Community for the fellowship 
no ERBCHBICT941569 "Multifractal Analysis of Spatio- Temporal Chaos". This work was 
partly supported by DEG within the GraduiertenkoUeg "Eeldtheoretische und numerische 
Methoden in der Elementarteilchen- und Statistischen Physik". Einally, we want to thank 
the Institute for Scientific Interchange (I.S.I) in Torino (Italy) for the hospitality offered 
us during July 1995 within a workshop organized by the EEC-Network "Complexity and 
Chaos". 



7 



REFERENCES 



[1] W.H. Press et al., Numerical Recipes (Cambridge Univ. Press, Cambridge 1986) 

[2] A. Kolmogorov, I. Petrovsky and N. Piskunov, Bull. Univ. Moscow^ Ser. Int. A 1 (1937) 

1; D.G. Aronson and H.F. Weinberger, Adv. Math. 30 (1978) 33. 
[3] M.C. Cross and P.H. Hohenberg, Phys. Rev. Mod. Phys. 65 (1993) 851. 
[4] P. Collet and J. P. Eckmann, Instabilities and Fronts in Extended Systems (Princeton 

University Press, Princeton, 1990) 
[5] A. Politi and A. Torcini, Europhys. Lett 28 (1994) 545. 
[6] A. Torcini, P. Grassberger and A. Politi, J. Phys. A 27 (1995) 4533. 
[7] L. Keefe, Phys. Lett. A 140 (1989) 317; B.I. Shraiman, A. Pumir, W. van Saarloos, 

P.C. Hohenberg, H. Chate and M. Holen, Physica D 57 (1992) 241; D. A. Egolf and 

H.S. Greenside, Nature 369 (1994) 129, Phys. Rev. Lett 74 (1995) 1751; T. Bohr, E. 

Bosch and W. van de Water, Nature 372 (1994) 48. 
[8] H. Chate, NonUnearity 7 (1994) 185. 
[9] P. Grassberger, Physica Scripta 40 (1985) 1033. 
[10] J. Crutchfield and K. Kaneko, Phys. Rev. Lett. 60 (1988) 2715; A. Politi, R. Livi, G.- 

L. Oppo and R. Kapral, Europhys. Lett 22 (1993) 571; A. Wacker, S. Bose and E. 

Scholl.Europhys. Lett. 31 (1995) 257. 
[11] J. Candy and W. Rozmus, J. Comp. Phys. 92 (1991) 230. 

[12] H. Yoshida, Phys. Lett. A, 150 (1990) 262; R.I. McLachlan and P.Atela, NonUnearity 
5 (1992) 541. 

[13] H. Erauenkron and P. Grassberger, Int. J. Mod. Phys. C 5 (1994) 37. 

[14] Higher order integration schemes, such as the fourth order method introduced by Candy 
and Rozmus [11], bring a noticeable reduction of errors for conservative Hamiltonian 
systems, but it turns out that such algorithms are not stable for dissipative systems such 
as the CLGE. This can be explained by the fact that in these higher order techniques 
also backward time steps are needed. 

[15] K. Kaneko, Prog, theor. Phys. 72 (1984) 980; I. Waller and R. Kapral, Phys. Rev. A 30 
(1984) 2047. 

[16] A. Politi and A. Torcini, Chaos 2 (1992) 293. 

[17] G.C. Paquette, L.Y. Chen, N. Goldenfeld and Y. Oono, Phys. Rev. Lett. 72 (1994) 76. 
[18] A. Torcini, P. Grassberger and H. Erauenkron, to be published. 



8 



FIGURES 



Fig.l: Real (solid line) and complex (dashed line) parts of the convolution kernel (j){x) for 
the CGLE (see eq.(13)), for ci = 3.5, Ax = 0.098 and r = 0.05. 

Fig. 2: Integration errors as a function of Nc for several integration time steps r. The errors 
are estimated by comparing with a standard pseudo-spectral code with a smaller time step 
t' = r/100. The values for Ci and Ax are the same as in fig.l, the other parameters are 
C3 = 0.9 and N = 2048. The errors shown are averages over several runs of duration T = 0.2 
each. 
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TABLES 



TABLE L Propagation speed and slope of the exponential profile of the propagating front for 
the CGLE with parameter ci = 3.5 at some values of C3. /^p and vp refer to a direct measurement 
through the algorithm here introduced, /^l and vp refer to the minimum of the function E(//) 
defined in Eq. (17). The direct measurement of the front propagation has been performed adopting 
the method here introduced with Nc = 110, A = 2048, r = 0.05, and Ax = 0.098. After a transient 
of 100,000 time steps, the front has been followed for a number of time steps ranging from 1,840,000 
to 2,100,000. The Lyapunov analysis necessary to estimate /^l and vp has been done employing a 
pseudo-spectral code with A = 2048, r = 0.01 and Ax = 0.098. After a transient of 100,000 time 
steps the dynamics of the system has been integrated for a number of steps ranging from 470,000 
to 1,290,000. In both cases periodic boundary conditions have been employed. 



C3 Vp ML 

(L9 ~ 0.16 2.09 ± 0.01 0.134 ± 0.006 2.13 ± 0.07 

0.8 ~ 0.12 1.34±0.01 0.110± 0.005 1.38±0.05 

0.7 ~ 0.10 0.95 ± 0.01 0.085 ± 0.005 0.96 ±0.04 

0.6 ~ 0.08 0.62 ± 0.01 0.075 ± 0.006 0.65 ±0.03 
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